
;;==========================================================================
;; Routine to make comparison plots for galactic positions derived
;;   from the NEAR/FAR determination of the code.  This analysis
;;   routine is a heavily-modified version of INVESTIGATE_GALPOS.pro

;; Define the COMMON BLOCK
COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
   lpstr,rho_hi,rho_h2,rho_star

;; Since we are using routines buried in BGPS_FFORE, compile now to
;; they are available.
FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_LOOPIE
RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER
 
;; Here's where you can change the SUBSET used (i.e. TEMP,
;; training set, etc.)  ALSO the SUFF is which rotation curve used
rotc   = 2b
subset = ''

do_Tfit = 0b

restore,'irdc_dist_model/bgps_rb3'+subset+'.sav',/ver
s = read_bgps_csv(/ver)
n = n_elements(rb3)


;; Define arrays for TEMP_DEPENDENT_VARIABLES
tau    = dblarr(n)
f_data = dblarr(n)
sig_f  = dblarr(n)

;; Get galactic params
defsysv, '!MW', exists = exists
IF NOT exists THEN galactic_params 

;; Using !MW parameters -- Goes into COMMON block!
R0 = !MW.R0
d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
corr = R0 / 8.5d3


message,'Setting the stage (re-used calculations)...',/inf
;; Do the calculations needed by every loop and place them into a
;; structure to save on computation time
lpstr = REPLICATE( CREATE_STRUCT('kdn',0.,'kdf',0.,$
                                 'r',dblarr(!MW.NBINS*5.),$
                                 'z',dblarr(!MW.NBINS*5.),$
                                 'rho_hi',dblarr(!MW.NBINS*5.),$
                                 'rho_h2',dblarr(!MW.NBINS*5.),$
                                 'rho_star',dblarr(!MW.NBINS*5.)), n)
;; Do distances for the particular curve
CASE rotc OF
   0: BEGIN                             ;; Reid et al. (2009)
      FOR i=0L, n-1 DO BEGIN
         
         lpstr[i].kdn = KDIST(rb3[i].l, rb3[i].b, rb3[i].vlsr, /NEAR)
         lpstr[i].kdf = KDIST(rb3[i].l, rb3[i].b, rb3[i].vlsr, /FAR) 
         
         omni_lbd2rz, rb3[i].l, rb3[i].b, d, r, z
         lpstr[i].R = r
         lpstr[i].Z = z
      ENDFOR
      suff = ''
   ENDCASE
   1: BEGIN                             ;; IAU standard
      FOR i=0L, n-1 DO BEGIN
         
         lpstr[i].kdn = KDIST_IAU(rb3[i].l, rb3[i].b, rb3[i].vlsr, /NEAR)
         lpstr[i].kdf = KDIST_IAU(rb3[i].l, rb3[i].b, rb3[i].vlsr, /FAR) 
         
         R0 = 8500.
         omni_lbd2rz, rb3[i].l, rb3[i].b, d, r, z, R0=R0
         lpstr[i].R = r
         lpstr[i].Z = z
      ENDFOR
      suff = '_iau'
   ENDCASE
   2: BEGIN                             ;; Clemens (1985)
      FOR i=0L, n-1 DO BEGIN
         
         lpstr[i].kdn = KDIST_CLEM(rb3[i].l, rb3[i].b, $
                                          rb3[i].vlsr, /NEAR)
         lpstr[i].kdf = KDIST_CLEM(rb3[i].l, rb3[i].b, $
                                          rb3[i].vlsr, /FAR) 
         
         R0 = 8500.
         omni_lbd2rz, rb3[i].l, rb3[i].b, d, r, z, R0=R0
         lpstr[i].R = r
         lpstr[i].Z = z
      ENDFOR
      suff = '_clem'
   ENDCASE
ENDCASE

;; Calculate the HI & H2 densities along the LOS for each source
;;   here, and place into the lpstr for later retrieval (speed-up).
FOR i=0L, n-1 DO BEGIN
   lpstr[i].rho_hi = hi_density(lpstr[i].R,lpstr[i].Z)
   lpstr[i].rho_h2 = h2_density(lpstr[i].R,lpstr[i].Z)
   ;; The 2400 = 2 * (H1 + H2) is the normalization
   lpstr[i].rho_star = 30.d * (exp( -(lpstr[i].R-R0)/2600.d - $
                                    abs(lpstr[i].Z)/300.d) + $
                               0.12 * exp( -(lpstr[i].R-R0)/3600.d - $
                                           abs(lpstr[i].Z)/900.d)) / (2400.d)
   
ENDFOR



;; Set up some of the variables needed throughout the program
chi = '!9' + string("143B) + '!X'
gr_rho = '!9' + string("162B) + '!X'
gr_tau = '!9' + string("164B) + '!X'
propto = '!9' + string("265B) + '!X'
;;"

;; There are currently 4 cases of the model
nmods     = 4
dfn       = strarr(nmods)
nearfar   = bytarr(nmods,n)
thisnf    = bytarr(n)
fitdist   = fltarr(n)
ffore_vec = fltarr(n)


;; Set up the plotting environment   
myps,'./irdc_dist_model/bgps_galpos_plot'+subset+suff+'.eps',$
         ct=0,xsize=12,ysize=6.5

cnames = ['YGB7','Lime Green','TG3','RED7']
colors = cgColor(cnames)
textcolor = cgColor('RYB8')
print,colors
thick=[6,5,4,3]

multiplot,[4,1],mtitle='Symsize '+propto+' IRDC Contrast',mtitoffset=-0.6

FOR Td = 10., 26, 5. DO BEGIN
   
   IF Td GT 21 THEN tstring = 'FIT' ELSE tstring = string(Td,format="(I0,'K')")
   
   message,'Dust temp = '+tstring,/inf
   
   ;; Put the temperature-dependent variables into the COMMON block.
   IF Td LE 21. THEN $
      TEMP_DEPENDENT_VARIABLES, Td, /BGPS
   
   IF Td EQ 10. THEN BEGIN
      ytit='Galactocentric Distance [kpc]'
      ytf = "(I0)"
   ENDIF ELSE BEGIN
      ytit=''
      ytf =''
   ENDELSE
   
   plotsym,0,2.0                ;,/fill
   plot,[-1.5,10.5],[-10,13.5],/nodata,xtit='Galactocentric Distance [kpc]',$
        ytit=ytit,ytickformat=ytf,/xst,/yst,tit='T!dd!n = '+tstring,$
        xtickformat="(I0)",charsize=0.75,/isotropic
   
   chmin = fltarr(nmods)
   ndof  = fltarr(nmods)
   tdp   = fltarr(nmods)
   
   FOR j=0, nmods-1 DO BEGIN
      
            message,'Case #'+string(j+1,format="(I0)"),/inf
      IF Td LE 21. THEN $
         dfn[j] = './irdc_dist_model/data/case'+$
                  string(j+1,format="(I0)")+'_irdc_Td'+$
                  string(Td,format="(I0,'K')")+subset+suff+'.sav' $
      ELSE $
         dfn[j] = './irdc_dist_model/data/case'+$
                  string(j+1,format="(I0)")+$
                  '_irdc_TdFIT'+subset+suff+'.sav'
      
      restore,dfn[j]
      chmin[j] = bestnorm
      ndof[j] = dof
      
      message,'Minimum Chisq = '+string(bestnorm,format="(F0.2)"),/inf
      
      IF Td GE 21 THEN do_Tfit = 1b ELSE do_Tfit = 0b
      
      residj = BGPS_FFORE_LOOPIE( res, NEARFAR=thisnf , FFORE_VEC=ffore_vec)
      
      ;; Sort out near vs. far vectors, and derive fitdist
      nearfar[j,*] = thisnf
      far = thisnf & near = ~ thisnf
      fitdist = near * lpstr.kdn + far * lpstr.kdf
      
      ;; Make x,y position for plotting
      x = fitdist/1.d3 * sin(rb3.l*!dtor)
      y = R0/1.d3 - fitdist/1.d3 * cos(rb3.l*!dtor)
      FOR kk = 0L, n-1 DO $
         plots,x[kk],y[kk],color=colors[j],psym=8,symsize=rb3[kk].c_meas
            
   ENDFOR
   
   ;; Solar Circle and Tangent Circles
   ells = findgen(1001)/1.e3*180.
   scl = findgen(1001)/1.e3*360.
   xells = (R0/1.d3 * cos(ells*!dtor)) *  sin(ells * !dtor)
   yells = - (R0/1.d3 * cos(ells*!dtor)) * cos(ells * !dtor) + R0/1.d3
   xsc = R0/1.d3 * cos(scl * !dtor)
   ysc = R0/1.d3 * sin(scl * !dtor)
   oplot,xells,yells,linestyle=1,color=cgColor('Olive'),thick=3
   oplot,xsc,ysc,linestyle=1,color=cgColor('Dark Slate Blue'),thick=3
   
   ;; l=30 line (through end of the Bar)
   plots,0,R0/1.d3
   plots,R0/1.d3*cos(30*!dtor),-R0/1.d3*sin(30*!dtor),/cont,linestyle=2,thick=3
   al_legend,position=[R0/1.d3*cos(30*!dtor),-R0/1.d3*sin(30*!dtor)],$
             box=0,charsize=0.75,['l = 30']
   
   ;; Sun and Galactic Center
   plots,0,0,psym=7,color=cgColor('Orchid'),symsize=1,thick=2
   plots,0,R0/1.d3,psym=7,color=cgColor('Sea Green'),symsize=1,thick=2
   
   plotsym,0,0.5,/fill
   al_legend,/top,/left,textcolor=textcolor,$
             ['Case '+string(indgen(nmods)+1,format="(I2)")+' ('+chi+$
              '!u2!dred!n = '+string(chmin/float(ndof),format="(F0.2)")+')'],$
             box=0,charsize=0.75,psym=8,colors=cnames
   
   multiplot
   
ENDFOR

multiplot,/reset
multiplot,/default
multiplot,/reset


device,/close_file

!p.font=-1
set_plot,'x'

END
